/*
 * Tcl interface for Geometric constraint using mbuild.
 */

#include "mbuild.h"
#include "tcl.h"
#include <ctime>
/*
 * Create Geo Constraint, current do not support save constraint data in a variable and only support ToroidalPore & GaussianConstraint.
 */
//TODO attach the external variables with a Tcl_Obj, and make these more generic
static int
CreateGeo_Cmd(ClientData cdata, Tcl_Interp *interp, int objc, Tcl_Obj *const objv[])
 {
     Tcl_Obj* geoVarObj=objv[1];
     Tcl_Obj* geoTypeObj =objv[2];
     ToroidalPore* pore;
     char result[256];
     char error[256];
     //Start ...
     double r0,k0,T,r_up,r_bottom,cutoff=10;
     int ncyc=3;
     if(objc <3) {printf("Wrong args, should be:\n CreateGeo name $type [args ..]\n");return TCL_ERROR;}
     /*test of using Tcl_obj internalRep*/
//     for(int i=0;i<objc;i++) {
//     printf("value: %s\n",objv[i]->bytes);
//     if (objv[i]->typePtr == NULL) printf("Unknown type\n");// typePtr is a static pointer
//     } 
//     printf("%d\n",Tcl_GetObjType("string"));
//    Tcl_ConvertToType(interp,geoVarObj,Tcl_GetObjType("string"));
//    Tcl_Obj* list_of_type=Tcl_NewObj();
//    Tcl_AppendAllObjTypes(interp,list_of_type);
//    printf("%s\n",Tcl_GetStringFromObj(list_of_type,NULL));
//    if (geoVarObj->typePtr != NULL) printf("type %s",geoVarObj->typePtr->name);
     /* get arguments and initialize*/
     // TODO : extendible geotypes?
     if(strcmp(geoTypeObj->bytes,"ToroidalPore")==0) {
         if(objc <6) {
             sprintf(error,"Wrong args, should be:\nCreateGeo name ToroidalPore $R0 $K0 $T [Rcut]\n");
             Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
             return TCL_ERROR;
         }
         Tcl_GetDoubleFromObj(interp,objv[3],&r0);
         Tcl_GetDoubleFromObj(interp,objv[4],&k0);
         Tcl_GetDoubleFromObj(interp,objv[5],&T);
         pore=new ToroidalPore();
         pore->setup(r0,k0,T);
         if(objc > 7) {
             Tcl_GetDoubleFromObj(interp,objv[6],&r_up);
             if(objc >8) Tcl_GetDoubleFromObj(interp,objv[7],&cutoff);
             if(objc >9) Tcl_GetIntFromObj(interp,objv[7],&ncyc);
             r_bottom=r_up;
             if(r_up >0) {
                 CylindricalHole hole=CylindricalHole(r_up,r_bottom,T+10*cutoff,-T/2-5*cutoff);
                  pore->has_hole(&hole);
                  }
         }
         pore->set_cutoff(cutoff);
         pore->set_integration_cyc(ncyc);
         Tcl_Obj* poreObj = Tcl_NewObj();
         Tcl_ObjSetVar2(interp,objv[1],NULL,poreObj,0);
         Tcl_ConvertToType(interp,geoVarObj,Tcl_GetObjType("bytecode"));
         (poreObj->internalRep).twoPtrValue.ptr1=pore;
         poreObj->bytes="ToroidalPore";
     } else if (strcmp(geoTypeObj->bytes,"EllipsePore")==0) {
         if(objc <6) {
             sprintf(error,"Wrong args, should be:\nCreateGeo name EllipsePore $Rin $Rout $T [Rcut]\n");
             Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
             return TCL_ERROR;
         }
         double r_in,r_out;
         Tcl_GetDoubleFromObj(interp,objv[3],&r_in);
         Tcl_GetDoubleFromObj(interp,objv[4],&r_out);
         Tcl_GetDoubleFromObj(interp,objv[5],&T);
         EllipsePore* pore2=new EllipsePore();
         pore2->setup(r0,k0,T);
         if(objc > 7) {
             Tcl_GetDoubleFromObj(interp,objv[6],&r_up);
             if(objc >8) Tcl_GetDoubleFromObj(interp,objv[7],&cutoff);
             if(objc >9) Tcl_GetIntFromObj(interp,objv[7],&ncyc);
             r_bottom=r_up;
             if(r_up >0) {
                  CylindricalHole hole=CylindricalHole(r_up,r_bottom,T+10*cutoff,-T/2-5*cutoff);
                  pore->has_hole(&hole);
                  }
         }
         pore2->set_cutoff(cutoff);
         pore2->set_integration_cyc(ncyc);
         Tcl_Obj* poreObj = Tcl_NewObj();
         Tcl_ObjSetVar2(interp,objv[1],NULL,poreObj,0);
         Tcl_ConvertToType(interp,geoVarObj,Tcl_GetObjType("bytecode"));
         (poreObj->internalRep).twoPtrValue.ptr1=pore2;
         poreObj->bytes="EllipsePore";
     } else {
         sprintf(error,"Unrecognized Geometric Type: %s!\n",geoTypeObj->bytes);
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
     }
     sprintf(result,"Created %s",geoTypeObj->bytes);
     Tcl_SetObjResult(interp, Tcl_NewStringObj(result, -1));
     return TCL_OK;
 }

/*
 * Set the Range of Geo Constraint, Inside this box, the constraint is added.
 */
static int
SetGeoRange_Cmd(ClientData cdata, Tcl_Interp *interp, int objc, Tcl_Obj *const objv[])
 {
    double xl,yl,zl,xh,yh,zh;
    char result[256];
    char error[256];
    //start ...
    if (objc < 2) {
           sprintf(error,"Wrong args, should be :\n setGeoRange name $low $high\n or \n setGeoRange name $xlow $xhigh $ylow $yhigh $zlow $zhigh\n");
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
    }
    Tcl_Obj* poreObj=Tcl_ObjGetVar2(interp,objv[1],NULL,0);
    if (poreObj == NULL ) {
        sprintf(error,"No variable %s exists in the current frame or namespace.",objv[1]->bytes);
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
    }
    if ((poreObj->internalRep).twoPtrValue.ptr1 == NULL ) {
         sprintf(error,"The variable is not in a corret type.");
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
    }
    Membrane* GeoShape = (Membrane*) (poreObj->internalRep).twoPtrValue.ptr1;
     if(objc <4 || objc % 2 != 0) {
         sprintf(error,"Wrong args, should be:\n setGeoRange name $low $high\n or \n setGeoRange name $xlow $xhigh $ylow $yhigh $zlow $zhigh\n");
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
     }
     if(objc == 4) {
         Tcl_GetDoubleFromObj(interp,objv[2],&xl);
         Tcl_GetDoubleFromObj(interp,objv[3],&xh);
         yl=zl=xl;
         yh=zh=xh;
     } else if (objc == 8) {
         Tcl_GetDoubleFromObj(interp,objv[2],&xl);
         Tcl_GetDoubleFromObj(interp,objv[3],&xh);
         Tcl_GetDoubleFromObj(interp,objv[4],&yl);
         Tcl_GetDoubleFromObj(interp,objv[5],&yh);
         Tcl_GetDoubleFromObj(interp,objv[6],&zl);
         Tcl_GetDoubleFromObj(interp,objv[7],&zh);
     } else {
         sprintf(error,"Wrong args, should be:\n setGeoRange name $low $high\n or \n setGeoRange name $xlow $xhigh $ylow $yhigh $zlow $zhigh\n");
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
     }
     if (strcmp(poreObj->bytes,"ToroidalPore")==0) {
         GeoShape->set_region(xl,xh,yl,yh,zl,zh);
     } else if (strcmp(poreObj->bytes,"EllipsePore")==0) {
         GeoShape->set_region(xl,xh,yl,yh,zl,zh,true);
     } else {
         GeoShape->set_region(xl,xh,yl,yh,zl,zh,true);
     }

     sprintf(result,"Box Range:\n x: %5f to %5f\n y: %5f to %5f\n z: %5f to %5f\n",xl,xh,yl,yh,zl,zh);
     Tcl_SetObjResult(interp, Tcl_NewStringObj(result, -1));
     return TCL_OK;
 }

/*
 * Generate Geo Constraint
 */
static int
GenGeoConstraint_Cmd(ClientData cdata, Tcl_Interp *interp, int objc, Tcl_Obj *const objv[])
 {
    char result[256];
    char error[256];
    int nx,ny,nz;
    double dx,dy,dz;
    nx=ny=nz=161;
    clock_t t1,t2;
    GridSpace* myConstraint;
    //start ...
     if(objc < 3) {
         printf("Wrong args, should be:\n genGeoConstraint name $type [args ...]\n");
         return TCL_ERROR;
     }
     // Load data from the variable given at first argument.
    Tcl_Obj* poreObj=Tcl_ObjGetVar2(interp,objv[1],NULL,0);
    if (poreObj == NULL ) {
        sprintf(error,"No variable %s exists in the current frame or namespace.",objv[1]->bytes);
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
    }
    if ((poreObj->internalRep).twoPtrValue.ptr1 == NULL ) {
         sprintf(error,"The variable is not in a corret type.");
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
    }
    Membrane* GeoShape= (Membrane*) (poreObj->internalRep).twoPtrValue.ptr1;
     
     /* get arguments and create constraint*/
     Tcl_Obj* consTypeObj=objv[2];
     // TODO : extendible consTypes?

    
     if(strcmp(consTypeObj->bytes,"Gaussian")==0) {
         double force,width,offset;
         if(objc >= 6) {
             Tcl_GetDoubleFromObj(interp,objv[3],&force);
             Tcl_GetDoubleFromObj(interp,objv[4],&width);
             Tcl_GetDoubleFromObj(interp,objv[5],&offset);
             if(objc >= 7) {
                          Tcl_GetIntFromObj(interp,objv[6],&nx);
                          ny=nz=nx;
                          dx=(GeoShape->xh-GeoShape->xl)/nx;
                          dy=(GeoShape->yh-GeoShape->yl)/ny;
                          dz=(GeoShape->zh-GeoShape->zl)/nz;
             }
         }   else {
             printf("Wrong args, should be:\n genGeoConstraint name Gaussian $forceconstant $width $offset [$ngrids]\n");
             return TCL_ERROR;
         }
         GaussianGeoConstraint* g=new GaussianGeoConstraint();
         try {
         g->set_forceconstant(force);
         g->set_width(width);
         g->set_offset(offset);
         } catch(ParameterError& error) {
             printf("Something wrong with the parameters(force,width,offset): %5f %5f %5f\n",force,width,offset);
             return TCL_ERROR;
         }
         GeoShape->set_distribution(g);
         myConstraint=new GridSpace();
         try {
         myConstraint->set_origin(GeoShape->xl,GeoShape->yl,GeoShape->zl);
         //Make the GridSpace size = (xh-xl)+dx,so that the margin at the high end could be used for linear interpolation.
         myConstraint->set_dimension(nx+1,ny+1,nz+1);
         myConstraint->set_grid_size(dx,dy,dz);
         } catch(ParameterError& error) {
             printf("Something wrong with the parameters:\n");
             printf("(x0,y0,z0) : %5f %5f %5f\n",GeoShape->xl,GeoShape->yl,GeoShape->zl);
             printf("(nx,ny,nz) : %5d %5d %5d\n",nx,ny,nz);
             printf("(dx,dy,dz) : %5f %5f %5f\n",dx,dy,dz);
             return TCL_ERROR;
         }
         myConstraint->initialize(force);
         t1=clock();
         GeoShape->fill(myConstraint,myConstraint);
         t2=clock();
         (poreObj->internalRep).twoPtrValue.ptr2=myConstraint;
     } else if(strcmp(consTypeObj->bytes,"Harmonic")==0) {
         double force,offset,cutoff;
         if(objc >= 6) {
             Tcl_GetDoubleFromObj(interp,objv[3],&force);
             Tcl_GetDoubleFromObj(interp,objv[4],&cutoff);
             Tcl_GetDoubleFromObj(interp,objv[5],&offset);
             if(objc >= 7) {
                          Tcl_GetIntFromObj(interp,objv[6],&nx);
                          ny=nz=nx;
                          dx=(GeoShape->xh-GeoShape->xl)/nx;
                          dy=(GeoShape->yh-GeoShape->yl)/ny;
                          dz=(GeoShape->zh-GeoShape->zl)/nz;
             }
         }   else {
             printf("Wrong args, should be:\n genGeoConstraint name Harmonic $forceconstant $cutoff $offset [$ngrids]\n");
             return TCL_ERROR;
         }
         HarmonicGeoConstraint* h=new HarmonicGeoConstraint();
         try {
            h->set_cutoff(cutoff);
            h->set_forceconstant(force);
            h->set_offset(offset);
         } catch(ParameterError& error) {
             printf("Something wrong with the parameters(force,offset,cutoff): %5f %5f %5f\n",force,offset,cutoff);
             return TCL_ERROR;
         }
         GeoShape->set_distribution(h);
         myConstraint=new GridSpace();
         try {
         myConstraint->set_origin(GeoShape->xl,GeoShape->yl,GeoShape->zl);
         //Make the GridSpace size = (xh-xl)+dx,so that the margin at the high end could be used for linear interpolation.
         myConstraint->set_dimension(nx+1,ny+1,nz+1);
         myConstraint->set_grid_size(dx,dy,dz);
         } catch(ParameterError& error) {
             printf("Something wrong with the parameters:\n");
             printf("(x0,y0,z0) : %5f %5f %5f\n",GeoShape->xl,GeoShape->yl,GeoShape->zl);
             printf("(nx,ny,nz) : %5d %5d %5d\n",nx,ny,nz);
             printf("(dx,dy,dz) : %5f %5f %5f\n",dx,dy,dz);
             return TCL_ERROR;
         }
         myConstraint->initialize(force);
         t1=clock();
         GeoShape->fill(myConstraint,myConstraint);
         t2=clock();
         (poreObj->internalRep).twoPtrValue.ptr2=myConstraint;
     } else {
         sprintf(error,"Unsupported constraint type!");
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
     }
     sprintf(result,"Created %s constraint in %dx%dx%d Grids for %5f seconds\n",consTypeObj->bytes,nx,ny,nz,double(t2-t1)/CLOCKS_PER_SEC);
     Tcl_SetObjResult(interp, Tcl_NewStringObj(result, -1));
     return TCL_OK;
 }
/*
 * Write Constraint Matrix into a Dx file
 */
static int
saveGeoConstraint_Cmd(ClientData cdata, Tcl_Interp *interp, int objc, Tcl_Obj *const objv[]) {
    char error[256];
    if(objc != 3) {
         printf("Wrong args, should be:\n saveGeoConstraint name $fname\n");
         return TCL_ERROR;
    }
    Tcl_Obj* poreObj=Tcl_ObjGetVar2(interp,objv[1],NULL,0);
    if (poreObj == NULL ) {
        sprintf(error,"No variable %s exists in the current frame or namespace.",objv[1]->bytes);
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
    }
    if ((poreObj->internalRep).twoPtrValue.ptr1 == NULL ) {
         sprintf(error,"The variable is not in a corret type.");
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
    }
    Membrane* GeoShape= (Membrane*) (poreObj->internalRep).twoPtrValue.ptr1;
    GridSpace* myConstraint= (GridSpace*) (poreObj->internalRep).twoPtrValue.ptr2;
    char* fname;
    int length;
    fname=Tcl_GetStringFromObj(objv[2],&length);
    try {
        DxFile::write_dx(fname,myConstraint);
    } catch(...) {
        printf("Cannot write to file %s !\n",fname);
        return TCL_ERROR;
    }
    Tcl_SetObjResult(interp, Tcl_NewStringObj(fname, -1));
    return TCL_OK;
}

/*
 * Get Constraint value
 *  getGeoConstraint geo $x $y $z v dvdx dvdy dvdz
 */
static int
GetGeoConstraint_Cmd(ClientData cdata, Tcl_Interp *interp, int objc, Tcl_Obj *const objv[])
 {
    double x,y,z,v,dvdx,dvdy,dvdz;
    Tcl_Obj* vObj;
    Tcl_Obj* dvdxObj;
    Tcl_Obj* dvdyObj;
    Tcl_Obj* dvdzObj;
    v=dvdx=dvdy=dvdz=0.0;
    char result[256];
    char error[256];
    //start ...
    
    Tcl_Obj* poreObj=Tcl_ObjGetVar2(interp,objv[1],NULL,0);
    if (poreObj == NULL ) {
        sprintf(error,"No variable %s exists in the current frame or namespace.",objv[1]->bytes);
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
    }
    if ((poreObj->internalRep).twoPtrValue.ptr1 == NULL || (poreObj->internalRep).twoPtrValue.ptr1 == NULL) {
         sprintf(error,"The variable is not in a corret type.");
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
    }
    Membrane* GeoShape= (Membrane*) (poreObj->internalRep).twoPtrValue.ptr1;
    GridSpace* myConstraint= (GridSpace*) (poreObj->internalRep).twoPtrValue.ptr2;
    //Debug information.
//    if (poreObj->typePtr != NULL) printf("%s type %s\n",poreObj->bytes,poreObj->typePtr->name);
//    printf("ptr1: %d ptr2: %d\n",GeoShape,myConstraint);
//    printf("poreObj value : %s\n",poreObj->bytes);
    
    /*must have name $x $y $z v */
     if(objc < 6) {
         printf("Wrong args, should be:\n getGeoConstraint name $x $y $z v [dvdx] [dvdy] [dvdz]\n");
         return TCL_ERROR;
     }
     /* get poistion*/
     Tcl_GetDoubleFromObj(interp,objv[2],&x);
     Tcl_GetDoubleFromObj(interp,objv[3],&y);
     Tcl_GetDoubleFromObj(interp,objv[4],&z);
     try {
         v=myConstraint->getDerivatives(x,y,z,dvdx,dvdy,dvdz,PEROIDIC);
     } catch(OutOfRange& e) {
         sprintf(error,"(%5f,%5f,%5f) is not in the range of the constraint box!",x,y,z);
         Tcl_SetObjResult(interp, Tcl_NewStringObj(error, -1));
         return TCL_ERROR;
     } catch(GridSpaceNotInitialized& error) {
         printf("GeoConstraint must be generated first!\n");
         return TCL_ERROR;
     }
     vObj=Tcl_ObjSetVar2(interp,objv[5],NULL,Tcl_NewDoubleObj(v),0);;
     /* if given variable name for returning dvdx,dvdy,dvdz then return the value*/
     if(objc>=9) {
         dvdxObj=Tcl_ObjSetVar2(interp,objv[6],NULL,Tcl_NewDoubleObj(dvdx),0);
         dvdyObj=Tcl_ObjSetVar2(interp,objv[7],NULL,Tcl_NewDoubleObj(dvdy),0);
         dvdzObj=Tcl_ObjSetVar2(interp,objv[8],NULL,Tcl_NewDoubleObj(dvdz),0);
     }
     sprintf(result,"%8f %8f %8f %8f",v,dvdx,dvdy,dvdz);
     Tcl_SetObjResult(interp, Tcl_NewStringObj(result, -1));
     return TCL_OK;
 }



/*
 * Initiate the package
 */
EXTERN int DLLEXPORT
 Geoconstraint_Init(Tcl_Interp *interp)
 {
    /* check if tcl version matched */
    /*
     if (Tcl_InitStubs(interp, TCL_VERSION, 0) == NULL) {
         return TCL_ERROR;
     }
     */
     /* changed this to check for an error  */
     if (Tcl_PkgProvide(interp, "GeoConstraint", "1.0") == TCL_ERROR) {
         return TCL_ERROR;
     }
     Tcl_CreateObjCommand(interp, "geoCreate", CreateGeo_Cmd, NULL, NULL);
     Tcl_CreateObjCommand(interp, "setGeoRange", SetGeoRange_Cmd, NULL, NULL);
     Tcl_CreateObjCommand(interp, "genGeoConstraint", GenGeoConstraint_Cmd, NULL, NULL);
     Tcl_CreateObjCommand(interp, "saveGeoConstraint", saveGeoConstraint_Cmd, NULL, NULL);
     Tcl_CreateObjCommand(interp, "getGeoConstraint", GetGeoConstraint_Cmd, NULL, NULL);
     return TCL_OK;
 }
